Analysis of monocyte data (SCZ vs CTR): complex model
Plotting the first QC results:
Density plot
Covariate correlation
#Covariate correlation
for (i in colnames(int.covs))
int.covs[[i]] <- as.numeric(int.covs[[i]])
MEs_lin(int.covs,int.covs)Next step: QC on vst-normalized expression values.
#Check influence of covariates on data variance after transformation
PCAplot(pca, "batch", PoV.df=PoV, colors=c('#2E7281','black','#50B8CF',"green"))
PCAplot(pca, "status", PoV.df=PoV)We next correct the data with the same model we use for the DESeq2 DEA (for plotting purposes post DEA).
The results are as follows: (1) PC-covariate correlation post-correction:
PCAplot(pca.corrected, "batch", PoV.df=PoV.corrected, colors=c('#2E7281','black','#50B8CF',"green"))
PCAplot(pca.corrected, "status", PoV.df=PoV.corrected)We then go into the differential gene-expression analysis:
out of 35316 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 812, 2.3% LFC < 0 (down) : 1257, 3.6% outliers [1] : 0, 0% low counts [2] : 0, 0% (mean count < 0) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results
As well as the number of differentially expressed genes at lfc </> -1/1 at padj < 0.05:
[1] 1807
EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
pCutoff = 0.1,
FCcutoff = 2,
labSize = 6,
xlim = c(-20,20),
ylim = c(-1,10),
legendPosition = 'right')#Plotting DEGs
ComplexHeatmap::pheatmap(batch.rem[rownames(dds) %in% top100.up,],
cluster_rows = T, cluster_cols = T, annotation_legend = T, show_colnames = F, show_rownames = F,
scale='row',
legend = T, treeheight_row = 0,
annotation_col=annot_df, annotation_colors = annot_colors, fontsize_row = 8)
ComplexHeatmap::pheatmap(batch.rem[rownames(dds) %in% top100.down,],
cluster_rows = T, cluster_cols = T, annotation_legend = T, show_colnames = F, show_rownames = F,
scale='row',
legend = T, treeheight_row = 0,
annotation_col=annot_df, annotation_colors = annot_colors, fontsize_row = 8)#DT::renderDT(data.frame(res.complex), "OH1.5A",scrollY=1000)
paste("Reminder: we have", length(upreg.genes), "upregulated and", length(downreg.genes), "downregulated genes.")[1] “Reminder: we have 787 upregulated and 1234 downregulated genes.”
dotplot(GO.up, split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale="free")
dotplot(GO.down, split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale="free")pheatmap::pheatmap(mic.panel, cluster_rows = T, cluster_cols = T, annotation_legend = F, show_colnames = F, show_rownames = T,
legend_breaks=c(-4,0,4,4), scale="row",
legend = F, treeheight_row = 0,
legend_labels=c("-4","0","4","Residual expression \n\n"),
annotation_col=annot_df, annotation_colors = annot_colors, fontsize_row = 8)
pheatmap::pheatmap(NFKB.panel, cluster_rows = T, cluster_cols = T, annotation_legend = F, show_colnames = F, show_rownames = F,
legend_breaks=c(-4,0,4,4), scale="row",
legend = F, treeheight_row = 0,
legend_labels=c("-4","0","4","Residual expression \n\n"),
annotation_col=annot_df, annotation_colors = annot_colors, fontsize_row = 8)
pheatmap::pheatmap(IFN.panel, cluster_rows = T, cluster_cols = T, annotation_legend = F, show_colnames = F, show_rownames = F,
legend_breaks=c(-4,0,4,4), scale="row",
legend = F, treeheight_row = 0,
legend_labels=c("-4","0","4","Residual expression \n\n"),
annotation_col=annot_df, annotation_colors = annot_colors, fontsize_row = 8)[1] “Microglia”
PCA_cov_cor_R(int.covs, mic.panel)
PCAplot(mic.pca, "status", PoV.df=NFKB.PoV, pc.1=5, pc.2=6)
paste("NFKB")[1] “NFKB”
PCA_cov_cor_R(int.covs, NFKB.panel)
PCAplot(NFKB.pca, "status", PoV.df=NFKB.PoV, pc.1=5, pc.2=6)
paste("IFNy")[1] “IFNy”